/************************************************************************
 * MechSys - Open Library for Mechanical Systems                        *
 * Copyright (C) 2012 Sergio Torres, Diego Fregolente                   *
 *                                                                      *
 * This program is free software: you can redistribute it and/or modify *
 * it under the terms of the GNU General Public License as published by *
 * the Free Software Foundation, either version 3 of the License, or    *
 * any later version.                                                   *
 *                                                                      *
 * This program is distributed in the hope that it will be useful,      *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the         *
 * GNU General Public License for more details.                         *
 *                                                                      *
 * You should have received a copy of the GNU General Public License    *
 * along with this program. If not, see <http://www.gnu.org/licenses/>  *
 ************************************************************************/


// Std lib
#include <math.h>

// MechSys
#include <mechsys/dem/domain.h>
#include <mechsys/util/fatal.h>
#include <mechsys/mesh/unstructured.h>
#include <mechsys/mesh/structured.h>

using std::cout;
using std::endl;
using std::ofstream;
using DEM::Domain;
using Util::PI;

struct UserData
{
    double                  A;       ///< Amplitude of the oscilation
    double                ome;       ///< Frequency of the oscilation
    double                 Rc;       ///< Radius of the cilynder for volume fraction calculation
    DEM::Particle         *Op;       ///< Pointer to the oscilating container
    Vec3_t                Xop;       ///< Initial position of the oscilating container
    DEM::Particle         *Op2;       ///< Pointer to the oscilating container
    Vec3_t                Xop2;       ///< Initial position of the oscilating container
    std::ofstream      oss_ss;       ///< file for the output data
};

void Setup (DEM::Domain & dom, void * UD) // this function will move the cylinder
{
    UserData & dat = (*static_cast<UserData *>(UD)); //Here the structure userdata assigned to the domain dom is recovered and can be used

    //Calculate the new position
    Vec3_t New_Pos = dat.Xop + dat.A*sin(dat.ome*dom.Time)*OrthoSys::e2; //dom.Time is the simulation time variable
    dat.Op->Position(New_Pos); //Putting the cylinder in positin New_Pos

    Vec3_t New_Pos2 = dat.Xop2 + dat.A*sin(dat.ome*dom.Time)*OrthoSys::e2; //dom.Time is the simulation time variable
    dat.Op2->Position(New_Pos2); //Putting the cylinder in positin New_Pos
    
    
    ///////////////////////////In here a viscous force is included//////////////////////////
    
    //double eta = 1.0;
    //double mu  = 1.5;
    //for (size_t i=0;i<dom.Particles.Size();i++) //This force goes trough every particle
    //{
        //DEM::Particle *P = dom.Particles[i];//Pointer to the ith particle
        //double Nrm       = norm(P->v);//The norm of the velocity vector
        //P->F += -eta*P->v*pow(Nrm,mu-1.0); // A force of the form -eta*velocity^mu
    //}
}

void Report (DEM::Domain & dom, void *UD)
{
    UserData & dat = (*static_cast<UserData *>(UD)); //Here the structure userdata assigned to the domain dom is recovered and can be used
    if (dom.idx_out==0) //opening the file for the first time
    {
        String fs;
        fs.Printf("%s_fabric.res",dom.FileKey.CStr());
        dat.oss_ss.open(fs.CStr());
        dat.oss_ss << Util::_10_6  <<  "Time" << Util::_8s << "Cn" << Util::_8s << "Phi" << std::endl;
    }
    if (!dom.Finished) 
    {

        //Calculating the volume fraction 
        Vec3_t Xmin,Xmax;
        dom.BoundingBox(Xmin,Xmax);
        double VolCon = (Xmax(2)-Xmin(2))*M_PI*dat.Rc*dat.Rc; // Xmax and Xmin are the limiting position of the free particles.
        double Phi = dom.Vs/VolCon;                   // Vs is a variable of domain containng to solid volume
//        printf("Phi e: %d\n" ,VolCon, Phi);
              
        //Calculating the coodination number Cn
        double Cn = 0;
        size_t n_inside = 0; // Number of particles not in contact with the boundaries
        for (size_t i=0; i<dom.Particles.Size(); i++)
        {
            if (!dom.Particles[i]->Bdry&&dom.Particles[i]->IsFree()&&(dom.Particles[i]->x(2)<Xmax(2)-1.5*dom.MaxDmax)) //The containers and particles in contact with them are ignored
            {
                Cn += dom.Particles[i]->Cn;
                n_inside++;
            }
        }
        Cn/=n_inside;
//        printf("Cn e: %d\n" ,Cn);
	
        dat.oss_ss << Util::_8s << dom.Time << Util::_8s << Cn << Util::_8s << Phi <<std::endl;
    }
    else dat.oss_ss.close();
}

int main(int argc, char **argv) try
{
    

    if (argc<2) throw new Fatal("This program must be called with one argument: the name of the data input file without the '.inp' suffix.\nExample:\t %s filekey\n",argv[0]);
    
    //number of threads
    size_t Nproc = 4; 
    if (argc==3) Nproc=atoi(argv[2]);
    String filekey  (argv[1]);
    String filename (filekey+".inp");
    if (!Util::FileExists(filename)) throw new Fatal("File <%s> not found",filename.CStr());
    ifstream infile(filename.CStr());
        

    //set the simulation domain ////////////////////////////////////////////////////////////////////////////
    //
    
    
    String savefile;    // File used on the box.cpp simulaiton without the .hdf5 suffix
    double verlet;      // Verlet distance for optimization
    size_t VOut;        // Decide is video should be render
    double Kn;          // Normal stiffness
    double Kt;          // Tangential stiffness
    double Gn;          // Normal dissipative coefficient
    double Gt;          // Tangential dissipative coefficient
    double Mu;          // Microscopic friction coefficient
    double dt;          // Time step
    double dtOut;       // Time step for output
    double To;          // Time span for oscillaiton
    double Tf;          // Final time for the test
    double Amp;         // Amplitude of oscillation
    double ome;         // 
    {
        infile >> savefile;     infile.ignore(200,'\n');
        infile >> verlet;       infile.ignore(200,'\n');
        infile >> VOut;         infile.ignore(200,'\n');
        infile >> Kn;           infile.ignore(200,'\n');
        infile >> Kt;           infile.ignore(200,'\n');
        infile >> Gn;           infile.ignore(200,'\n');
        infile >> Gt;           infile.ignore(200,'\n');
        infile >> Mu;           infile.ignore(200,'\n');
        infile >> dt;           infile.ignore(200,'\n');
        infile >> dtOut;        infile.ignore(200,'\n');
        infile >> To;           infile.ignore(200,'\n');
        infile >> Tf;           infile.ignore(200,'\n');
        infile >> Amp;          infile.ignore(200,'\n');
        infile >> ome;          infile.ignore(200,'\n');
    }

    // Initializing the User data structure for the oscillating cylinder
    UserData dat;
    DEM::Domain dom(&dat); // This constructor will associate the UserData dat with the domain dom
    
    // Camera position
    double cam_x=0.0, cam_y=180.0, cam_z=0.0;
    dom.CamPos = cam_x, cam_y, cam_z;
    dom.Alpha=verlet;

    //Loading the sample prepared by box.cpp
    dom.Load(savefile.CStr());

    //Deleting all the obstacles since they are not neccesary for this part
    Array<int> DeleteTags(11);
    DeleteTags = -3,-4,-5,-6,-7,-8,-9,-10,-11,-12,-15;   // These tags are of all the obstacles
    //Array<int> DeleteTags(6);
    //DeleteTags = -2,-3,-4,-5,-6,-7;
    dom.DelParticles(DeleteTags);

    // Centering the free particles
    dom.Center(OrthoSys::O);

    //Creating the cylinder and the plane
    Vec3_t Xmin,Xmax; //The maximun and minimun positions of the rectangular box containing the particles.
    dom.BoundingBox(Xmin,Xmax);
    double Rc = 0.5*std::max(Xmax(0)-Xmin(0),Xmax(1)-Xmin(1))+dom.Particles[0]->Props.R; 
    //double Rc = sqrt(pow(Xmax(0)-Xmin(0),2.0) + pow(Xmax(1)-Xmin(1),2.0))+dom.Particles[0]->Props.R; 
    dat.Rc    = Rc;  //Storing the radius of the cylinder in the UserData structure for the volume fraction calculation
    Vec3_t Xc1(0.0,0.0,Xmin(2)-dom.Particles[0]->Props.R);
    Vec3_t Xc2(0.0,0.0,6.0*Xmax(2));

    dom.AddCylinder (/*int Tag*/-2, /*Vec3_t const & X0*/Xc1, /*double R0*/Rc, /*Vec3_t const & X1*/Xc2, /*double R1*/Rc, /*double R*/dom.Particles[0]->Props.R, /*double rho*/10.0*dom.Particles[0]->Props.rho);
    dom.AddPlane    (/*tag*/-3,/*position*/Xc1,/*spheroradius*/dom.Particles[0]->Props.R,/*Lx*/100,/*Ly*/100,/*rho*/10.0);
    
    //Fixing the obstacle particles
    dom.GetParticle(-2 ,true)->FixVeloc();
    dom.GetParticle(-3 ,true)->FixVeloc();


    // Setting the parameters
    Dict B;
    B.Set(-1 ,"Kn Kt Gn Gt Mu",Kn,Kt,Gn,Gt,Mu);
    B.Set(-2 ,"Kn Kt Gn Gt Mu",Kn,Kt,Gn,Gt,Mu);
    B.Set(-3 ,"Kn Kt Gn Gt Mu",Kn,Kt,Gn,Gt,Mu);
    dom.SetProps(B);

    //The cylinder (tag -2) will oscillate
    dat.ome = 2*M_PI*ome/To; //The frequency will be stored in the UserData strcture
    dat.A   = Amp; //The amplitude of the force will be stored in the UserData strcture
    dat.Op  = dom.GetParticle(-2,true); //The pointer to the cylinder is stored into the structure
    dat.Xop = dat.Op->x; //The initial center of mass positin of the cylinder will be stored for future reference
    //dat.Op->vxf = false; //Allow the particle to move along the x axis.
    
    //The plane (tag -3) will oscillate
    dat.ome = 2*M_PI*ome/To; //The frequency will be stored in the UserData strcture
    dat.A   = Amp; //The amplitude of the force will be stored in the UserData strcture
    dat.Op2  = dom.GetParticle(-3,true); //The pointer to the cylinder is stored into the structure
    dat.Xop2 = dat.Op2->x; //The initial center of mass positin of the cylinder will be stored for future reference
    //dat.Op->vxf = false; //Allow the particle to move along the x axis.

    
    
    //Initialize the gravity on the particles
    for (size_t i=0;i<dom.Particles.Size();i++)
    {
        dom.Particles[i]->Ff = dom.Particles[i]->Props.m*Vec3_t(0.0,0.0,-9.8);
    }

    String fkey_a(filekey+"_o"); //The files with the suffix 'o' are for the oscillation part
    dom.Solve     (/*tf*/To, /*dt*/dt, /*dtOut*/dtOut, &Setup, NULL, fkey_a.CStr(),VOut,Nproc);
    dom.Save(fkey_a.CStr());

    // After the oscillation the cylinder is no longer shaked hence the funciton setup is not called
    String fkey_b(filekey+"_r"); //The files with the suffix 'r' are for the relaxation part
    dom.Solve     (/*tf*/To+Tf, /*dt*/dt, /*dtOut*/dtOut, NULL, &Report, fkey_b.CStr(),VOut,Nproc);
    dom.Save(fkey_b.CStr());

    return 0;
}
MECHSYS_CATCH
